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FOREWORD 


This document contains a collection of eight papers presented at the AIAA/ASME/ASCE/ 
AHS/ASC 30th Structures, Structural Dynamics and Materials Conference held in Mobile, 
Alabama, April 3-5, 1989. The conference had a total of 278 papers, including 249 full-length 
papers, 17 papers in special and panel sessions, and 12 short presentations in the work-in- 
progress sessions. Seven of the papers appearing in this document were presented in the two 
work-in-progress sessions and the eighth paper was presented in one of the panel sessions. Most 
of the full length papers are contained in the conference proceedings published by AIAA. 

The fields covered by the conference are rapidly changing, and if new results and antici- 
pated future directions are to have maximum impact and use, it is imperative that they reach 
workers in the field as soon as possible. This consideration led to the decision to publish these 
proceedings prior to the conference. Special thanks go to the Research Information and Applica- 
tions Division at NASA Langley Research Center for their cooperation in publishing this 
volume. 

The use of trademarks or manufacturers’ names does not constitute endorsement, either 
expressed or implied, by the National Aeronautics and Space Administration. 


William F. Hunter 
Ahmed K. Noor 

Compilers 
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MODEL DESCRIPTION 


Delamination characterization of a plate with a hole/inclusion is presented. 
A closed form solution is developed to obtain stresses on the boundary of 
hole/inclusion in the plate. Once the location with highest tangential stress 
is identified, FEM analysis of a laminate under tensile loading is considered. 
The models for the closed form and FEM are given in Figure 1 A & B. The 
examples are for [±35/0/90] g laminate of AS4/3501-6. Tensile, Biaxial and Shear 
loads are considered. 


Y 



Figure 1 



PROCEDURE 


The analysis procedure is summerized in the following figure. The 
combination of closed-form (CFH - closed-form, hole [1] and CHI - closed-form, 
inclusion [2]) and PEM [3] solutions provide efficiency and economy in the 
interpretation of results. 


Note: CLT is a Classical Laminate Plate Theory Algorithm 



Closed 

Form 

Solution 


4 - 


Strain Energy Release Rates 


Figure 2. 
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INPLANE STRAIN/ STRESS RESULTS 


{±35/0/90] laminate is subjected to 5000 psi axial stress. The response of 
this laminate with a hole and with inclusions are tabulated below. The 
tangential stress is obtained from the closed-form solution and the strain is 
evaluated from CLT. 



E inclusion/ 
E plate 

V 

a applied 

t 

o 

hole 

(w/o inclusion) 

0 

3.280 

1710ye 

with soft 
inclusion 
(epoxy) 

.169 

1.865 

972ye 

with rigid 

inclusion 

(steel) 

3.128 

.724 

377ye 


Figure 3 
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TANGENTIAL STRESS AROUND THE HOLE BOUNDARY 


The tangential stress field for a laminate with a hole subjected to uniaxial 
tension and shear is presented in Figure 4. These stresses are displayed as a 
function of 9 around the boundary of the hole. As can be observed, the highest 
tangential stress is at 9 = 90° for the tensile load and 9 = 50.5° for the shear 
load. 



Figure 4 
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INTERLAMINAR NORMAL STRESS VARIATION 


The closed-form reslults of Figure 4 are converted to applied tensile 
strain, e * 1710pe to be used in the FEM calculations. The FEM analysis 
models th8 X laminate cross-section through-the-thickness and evaluates the 
interlaminar normal and shear stresses. Figure 5 presents the variation of the 
interlaminar normal stress along the y axis for each interface. Note that high 
stresses are at the free edge, between the 90/90 laminae and the 0/90 laminae. 


LEGEND 

INTERFACE 90/90 

- - INTERFACE 0/90 
- - INTERFACE 0/— 35 
INTERFACE -35/35 


^ 1 000 • 


-1000 


0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.1 

Y AXIS, in 


Figure 5 




INTERLAMINAR SHEAR STRESSES 

Figure 6 presents the variation of the interlaminar shear stresses, <j zx and 
a , along the Y-axis for each interface. It should be observed that the 
hf^hest a is at the free edge between the -35/35 laminae and that a rapidly 
vanishes fnside the laminate. Note that the highest a is just insiie the 
laminate near the free edge between the 0/-35 laminae. 2 ^ 
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INTERLAMINAR STRESSES THROUGH THE THICKNESS 


Figure 7 presents a comparison of the relative magnitudes of a , a 

at the free edge through the thickness. As can be observed, <x is much smaller 
than a and a . a is relatively large at the midplane, 90/90, and at the 
first interface, 0/90. The shear stress is highest at the third interface, 
-35/35. The fourth interface is the top surface of the laminate and has no 
interlaminar stresses. 

The strain energy release rates for this laminate are evaluated by modeling 
a crack of a length equal to the thickness of 8 plies at the 0/90 interface. 

The resulting ratios are G T / rTf)T = .94096 and G TT / rTnT * .05827, indicating that 
Mode I behavior is dominant. 



Figure 7 






DEFORMATION OF HOLE/INCLUSION BOUNDARY 


The effects of inclusion moduli on deformation are displayed in Figure 8. 
As expected a soft inclusion results in larger displacements. The [+35/0/90] 
laminate is loaded with uniaxial tensile stress of 50 ksi. The inclusion 
diameter is .5 inches. 



Figure 8 
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INTEGRALLY MACHINED WAFFLE PANELS 


Waffle panels are often used on fuselage shell structures such as that of the 
Space Shuttle. There are a number of advantages to using waffle paneling. The 
waffle panel design is an efficient design for carrying biaxial, in-plane and 
shear loads. The geometry of the waffle panel pockets can be adjusted to account 
for variability of loading in the principle directions. The stiffeners of the 
waffle panel readily provide attachment support for secondary systems, an impor- 
tant consideration for aerospace structures. The waffle panel also represents the 
simplest integrally machined panel design to manufacture. Figure 1 shows a typi- 
cal waffle panel design. The integrally machined waffle panel is constructed with 
a fillet radius between the vertical stiffener and the plate. This fillet radius 
provides additional stability and load carrying capability to the structure. 
Present analysis techniques for waffle type structures include classical theory of 
plate buckling, introduced here as the Rockwell-developed program WAFFLE, and the 
NASA-developed stress analysis program — Panel Analysis and Sizing Code (PASCO). 
This paper discusses the application of the PASCO program in conjunction with the 
WAFFLE program to account for both the fillet radius and the presence of 
stiffeners in both directions. The results of the tests are used to verify that 
these adjustments are valid and necessary if accurate analysis of the waffle panel 
is to be achieved. 



Figure 1 
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THE WAFFLE PROGRAM 


The WAFFLE program was developed at Rockwell for application on waffle 
panels. It is based on NACA Technical Notes (refs. 1, 2, 3, 4, and 5). The 
critical stresses are determined on the basis of the principle that during 
buckling, the elastic strain energy stored in a structure is equal to the work 
done by the applied loads. The deflection function is expressed exactly by means 
of a two-dimensional infinite Fourier series. The Raleigh-Ritz method is applied 
to obtain an infinite set of homogeneous linear equations needed to solve for the 
Fourier coefficients. The solutions of these equations that give Fourier coeffi- 
cients not all equal to zero exist only for those combinations of shear and direct 
stress for which the buckled plate is in neutral equilibrium. The matrix itera- 
tion method was then applied to solve for the 10 most important equations that 
best satisfy the loading condition. Input load includes in-plane and moment loads 
in both directions, shear loading, and lateral pressure. Thermal stress caused by 
temperature gradients is also calculated. The effect of the combined loadings is 
reflected in the general instability margin of safety. Local instability of the 
waffle panel is also taken into account by computing the crippling and buckling 
allowables of the stiffener as well as pocket skin buckling allowable (ref. 6). 

It should be noted that although it is the same plate buckling theory that is used 
to obtain the buckling allowable of the whole panel as well as the buckling allow- 
able of the pocket skin, the bending stiffness formulation of the former reflects 
the orthotropic effects presented by the biaxial stiffeners. WAFFLE therefore 
represents a comprehensive stress analysis of the waffle panel. Figure 2 shows a 
sample output of the program. 


COMPRESSION TEST 

PANEL GEOMETRY AND PROPERTIES AT 70. DEG. TEMPERATURE: 


A =20.00 

TS =.0660 
TEX =.0921 
RAD =1875 

B =18.00 
H =.800 
TEY =.0856 
FCY =57037. 

BSX =4.00 
TWX= .086 
IX =.003530 
E= 10.9E + 06 

BSY =3.00 
TWY= .086 
IY =.002775 

APPLIED LOADS: 

( + IS COMPRESSION AND - IS TENSION) 

NX =2000. 
TEMP= 70.0 

NY = 0. 

TDX = .00 

NXY = 0. 

TDY = .00 

PRESS = .00 

POCKET STABILITY: (LBS/IN.) 



NXP = 1 434. 

NYP =0. 

NXYP = 0. 

M.S. = -.03 

STIFFENER STABILITY: (PSD 



FCSX = 21 722. 
FCSY = 0. 

FCWX = 21722. 
FCWY = 0. 


M.S. = 1 .48 
(NO LOADS) 

GENERAL STABILITY: 



GRX = .9010 
RTX = .0000 

GRXY.0000 
RTV .0000 


M.S. = s.11 


Figure 2 
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THE PASCO PROGRAM 


The PASCO program was developed by M.S. Anderson and W.J. Stroud of NASA. It 
was designed for analyzing and sizing uniaxially stiffened panels. Buckling and 
vibration analyses are carried out with a linked plate analysis computer code 
denoted VIPASA, which is incorporated into PASCO (refs. 7, 8, and 9). Typical 
loading conditions for an uniaxially stiffened panel include longitudinal and 
transverse loads, shear load, bending moment, lateral pressure, temperature, and a 
bow-type imperfection (fig. 3). PASCO ordinarily models a cross section by 
assuming there is a repetition of substructures. The substructure is composed of 
an arbitrary assemblage of thin, flat, rectangular plate elements that are con- 
nected together along their longitudinal edges. The loads on each plate element 
are calculated under the assumption of uniform longitudinal strain. In addition, 
transverse load is assumed to be carried by the skin elements. 
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PASCO MODELS 


PASCO's capability is limited to uniaxially stiffened panels. Biaxially 
stiffened panels such as waffle panels, however, can be analyzed by assuming that 
the waffle is made up of a series of uniaxially stiffened panels whose interface 
lines are the location of the transverse stiffeners. Boundary conditions along 
the longitudinal edges can be simply supported or clamped. Boundary conditions 
can not be prescribed on the ends (transverse edges) of the panel. PASCO, how- 
ever, implicitly assumes that the end lines remain straight after buckling. Thus, 
in effect, PASCO imposes a simply support condition on the ends. Since the real 
waffle panel behaves neither as simply supported or clamped edge condition at the 
interface lines, we would expect a conservative result. A more accurate model 
would apply a bending moment at the interface. This bending moment can be deter- 
mined from WAFFLE program. Additional loading because of the thermal gradient 
could indirectly be determined by WAFFLE and be input into PASCO. Two PASCO 
models were run and the results were compared with the actual test result. Panel 
geometry, loading, and boundary conditions similar to the actual test were made. 

In one of the PASCO models, five small elements stacked on top of each other sim- 
ulated the fillet radii. The total area of these elements added up to the total 
area of the real fillet. Plots of the two models are shown in figure 4. 




Figure 4 
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EXPERIMENTAL 


Two full scale tests were run on the integrally machined waffle panels 
(ref. 10). For both tests, a series of uniaxial strain gages were used to deter- 
mine the stress in the expected area of failure. The gages were placed in the 
pockets and, where applicable, on the base, sides, and tip of the stiffeners. The 
first test was a study of the pocket buckling allowable for the waffle skin 
panel. In this test the stiffener was machined so that the skin was critical. 

The second test was a study of the stiffener buckling allowable. In this test, 
the stiffener was undercut to make the stiffener critical. The test panels were 
19 by 20 inches and 14 by 20 inches, respectively. They were installed on a 
Tinius Olsen Compression Test Machine, which can deliver up to 440,000 pounds of 
load. The load was applied in the longitudinal direction. The panels were placed 
vertically on the test machine* A metal slab sandwiched between the panel ends 
and the test bed uniformly distributed the compressive load. In order to simulate 
simply support condition, I-beams were firmly placed, but not clamped, along the 
longitudinal edges. In both tests the load was increased incrementally by 5,000 
pounds to the point of failure with stress readings taken at the end of each 
interval. The stress readings were tabulated with the corresponding applied 
loads. Geometry dimensions used in the two tests are shown in figure 5. A picture 
of the test setup is shown in figure 6. 



GEOMETRY 

TEST 1 

TEST 2 

B 

3.00 IN. 

2.24 IN. 

H 

0.80 IN. 

1 .13 IN. 

Ts 

0.066 IN. 

0.070 IN. 

Tw 

0.086 IN. 

0.058 IN. 

R 

0.19 IN. 

0.19 IN. 

L 

4.00 IN. 

3.38 IN. 

a x b 

1 

19 x 20 IN. 

14 x 20 IN. 


Figure 5 
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TEST SETUP 


ORIGfNAL PAGET 

BLACK AND WHITE PHOTOGRAPH 



Figure 6 
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RESULTS 


Plots of stress versus longitudinal compressive load are shown in figure 7. 
The relationship between the applied load and the induced stress is, as expected, 
very close to linear. Stress readings were taken from the appropriate strain 
gages mounted in the buckled areas (pocket skin in the first test and the 
stiffener in the second test). The results of the tests as well as the results 
predicted by WAFFLE and PASCO are shown in the following chart. The actual test 
result in each test is the average of the two strain gages. 



Test 1 

Buckling Half-Wavelength 

Test 2 

Buckling Half-Wavelength 

Method 

X = L/2 

X = L 

X = L/2 

X = L 

WAFFLE 

N/A 

21,121 psi 

N/A 

15,259 psi 

PASCO without fillet 

24,492 psi 

23,127 psi 

26,316 psi 

31,205 psi 

PASCO with fillet 

28,805 psi 

33,319 psi 

38,067 psi 

43,473 psi 

Actual test 

30,000 psi 

N/A 

39,750 psi 

N/A 



FIGURE 7. RESULTS 



FIGURE 7. RESULTS (CONTI 


Figure 7 
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INTRODUCTION 


DYSCO is a versatile, general purpose dynamic analysis program which assembles 
equations and solves dynamics problems. The "executive" manages a library of 
technology modules which contain routines that compute the matrix coefficients of 
the second order ordinary differential equations of the components. The executive 
performs the coupling of the equations of the components and manages the solution of 
the coupled equations. 

Any new component representation may be added to the library if, given the 
state vector, one can write a FORTRAN program to compute M, C, K, F. The problem 
described in this report demonstrates the generality of this statement. 


o DYSCO - MODELS AND SOLVES DYNAMIC SYSTEMS 

o EXECUTIVE MANAGES A LIBRARY OF TECHNOLOGY MODULES 

0 TECHNOLOGY MODULES DEFINE SECOND ORDER ODE <M, C, K, F) 

0 EXECUTIVE COUPLES COMPONENTS AND MANAGES SOLUTIONS 

0 NEW MODULES MAY BE ADDED TO LIBRARY IF: 

GIVEN THE STATE OF THE COMPONENT. A 
FORTRAN PROGRAM CAN COMPUTE M. C. K. F 

Figure 1 
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MISSILE - LAUNCH TUBE PROBLEM 


Consider an elastic body (missile) moving through an elastic tube. The inside 
of the tube is populated with snubbers having spring and damper characteristics with 
gaps between the missile and the snubbers. 

The complete problem also must include aerodynamic effects, control devices and 
algorithms, and propulsive force history. All of these are readily treated in 
DYSCO. 

The emphasis of this presentation will be the treatment of time varying 
constraints. 



CHARACTERISTICS 

MISSILE AND TUBE ARE ELASTIC BODIES 
GAPS EXIST BETWEEN MISSILE AND SNUBBERS 
AERODYNAMIC EFFECTS, CONTROLS, PROPULSIVE FORCE MUST BE 
CONSIDERED 

BASIC PROBLEM IS TIME VARYING CONSTRAINTS 


Figure 2 
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NEW DYSCO COMPONENT, CMS1 


The component illustrated in Figure 3a was added to the technology library. 
The "MS" stands for "moving structure." Note that there are no dynamic charac- 
teristics associated with the degrees of freedom in this component. The names of 
the DOF Z, , ... Z w are arbitrary and supplied by the user. When these correspond to 

t h 0 S G ^ ^ 

of another component, DYSCO automatically performs the dynamic coupling. The 
same is true of the DOF Z^ and Z^. 

Figure 3b illustrates a particular state of the system and illustrates how the 
force vector is coupled (M, C, K are treated as null). In this illustration, forces 
act only on DOF Z4, Z5, and ZU. All other forces are zero. 



Figure 3a 




Figure 3b 
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MISSILE-LAUNCH TUBE MODEL 


The figure illustrates the DYSCO model used to represent this problem. The 
missile and the tube are represented by modal representations (CFM2), six snubbers 
as previously described (CMS1) are used, each identified by the tube station at 
which it is located. 



****************************** 
BASIC MISSILE IN TUBE 

MODEL MT 

INDEX 

COMP NO. 

DATA SET 

FORCE 

1 

CFM2 1 

TUBE 

NONE 

2 

CFM2 2 

MSSL 

NONE 

3 

CMS1 

STA10 

NONE 

4 

CMS1 

STA20 

NONE 

5 

CMS1 

STA30 

NONE 

6 

CMS1 

STA40 

NONE 

7 

CMS1 

STA50 

NONE 

0 

CMS1 

STA60 

NONE 

A A A*A ****************************** 

************** *i 


********** ****************** 


DATA SET 


Figure 4 
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TUBE DATA SET 


This set of data defines the tube used in the sample analysis. Note the DOF 
names shown are also used on the snubber springs and dampers and are thus auto- 
matically attached to the tube at the proper station. 


o 


ZL10 ZL20 ZL30 


10" 20* 30" 


ZL40 ZLSO ZL60 


40" 50" 60" 


1 elastic node, 5 h z 


»*«,*•*«**»******««««***»«.«. TUBE / CFM2 

MISSILE LAUNCHER - NO X DOF. COMP 1. 


INPUT FOR STRUCTURAL COMPONENT CFM2 . MODAL FUSELAGE 


1 RBM 

2 NHOO 

3 NS 

4 X 


5 VC1 

6 Zi 


7 ZP1 


8 LC1 

9 TCI 

10 NR 

11 N I 

12 CIDFL 


13 XSTA 


14 AF 


15 AL 


16 MMS 

17 MO 
IB FREQ 


- RIGID BODY MODES - NO 

- NO. OF ELASTIC MOOES • l 

- NO. FUSELAGE STAS - 6 

- (REAL) INPUT STATION VALUES 

2 . 00000£+01 - 1 . 0000QE+01 0. OOOOOE+OO 1.000G0E+01 

2.00000E+01 3 . OOOOOE+Ol 

- H0DE1 VERTICAL COMP - YES 


- (REAL) MODE 1 VERTICAL DISP 

2 . OOOOOE -01 -4 . OOOOOE'Ol -6.00000E-01 -4.00000E-01 

2 . OQOOOE-Ol 1 .OOOOOE+OO 

- (REAL) MODE 1 VERTICAL SLOPE 

0 . OOOOOE+OO 0 . OOOOOE+OO 0. OOOOOE+OO 0. OOOOOE+OO 
0 . OOOOOE+OG 0 . OOOOOE+OO 


- MODE 1 LATERAL COMP - NO 

- MODE l TORSION COMP - NO 

- NO. OF ROTORS - 0 

- NO. OTHER IMPLCT DOF * 6 


- (DOF) IMPLICIT DOF NAMES 

ZL 10 0 ZL20 0 ZL30 0 ZL40 0 ZLSO 0 

ZL60 0 


- (REAL) STAS FOR IMPLCT DOF 
- 2 . OOOOOE+Ol -1. OOOOOE+Ol O.QGOOOE+OI 1. OOOOOE+Ol 

2 . OOOOOE+Ol 3 . OOOOOE+Ol 

- (REAL) FWO ANGLE FROM VERT 

0 . OOOOOE+OO 0 . OOOOOE+OO 0 . OOOOOE+OO 0. OOOOOE+OO 

0. OOOOOE+OO 0. OOOOOE+OO 

- (REAL) LAT ANGLE FROM VERT 

0. OOOOOE+OO 0 . OQOQOE+OQ 0. OOOOOE+OO 0. OOOOOE+OO 

0. OOOOOE+OO 0. OOOOOE+OO 


- MODAL MASS (SLUGS) • 2. OOOOOE+OO 

- MODAL DAMPING (PCI) - 0. OOOOOE+OO 

- MODAL FREQUENCY (HZ) - 5. OOOOOE+OO 


Figure 5 
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MISSILE DATA SET 


This set of data defines the missile used in the sample analysis. Note the DOF 
names which are used in the new component. 


ZO Z5 ZIO Z15 Z20 



0 5" 10" 15“ 20“ 



ISSUE - 

X OOF. COMP 2. 

HSSL 

/CFM2 ** 

NPUT FOR 

w BB JB ■ B B M Jo m ^ 

STRUCTURAL COMPONENT 
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STA 10 SNUBBER DATA SET 


Note the relationship between the DOF and those of the missile and the tube on 
previous figures. 
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STA 60 SNUBBER DATA SET 


Note DOF names on this data set and how they relate to missile and tube data. 



ZL60 


zo 


zs 




XCG 2 
XS - 50- 


i 



ZL60 


TUBE SNUBBER. STA 60 


STA60/CMS1 


INPUT FOR COHPONENT CMS1. H0VIN6 CONSTRAINTS 


1 

iU 

- 

UPPER BASE OOF 

m 

YES 

2 

COFLU 

- 

UPPER BASE OOF NAHE 

m 

ZL60 0 

3 

KU 

- 

UPPER SPRING COEFF 

» 

2 . OO0OOE+O2 

4 

CU 

- 

UPPER DAMPER COEFF 

m 

0. OOOOOE+0O 

5 

GU 

- 

UPPER GAP 

m 

5 . OOOOOE-Ol 

6 

1L 

- 

LOWER BASE DOF 

m 

YES 

7 

COFLL 


LOWER BASE DOF NAME 

M 

ZL60 0 

8 

KL 

- 

LOWER SPRING COEFF 

m 

2 . 00000E+02 

9 

CL 

- 

LOWER DAMPER COEFF 

m 

0. OOOOOE+OO 

10 

GL 

- 

LOWER GAP 

* 

5. OOOOOE-Ol 

11 

COFLX 

- 

HORIZONTAL DOF NAHE 

■ 

XCG 2000 

12 

XS 

- 

BASE DISTANCE 

m 

5. OOOOOE+Ol 

13 

NCOFZ 

- 

NO. OF VERTICAL DOF 

m 

s 

14 

CDFLZ 

ZO 

(DOF) VERTICAL DOF NAMES 
0 ZS 0 Z10 0 Z15 0 Z20 

15 

XZ 

- 

(REAL) VERT DOF STATIONS 


-1 . 00000E+01 -S.OOOOOE+OO O.OOOOOE+OO S.OOOOOE+OO 
1 . 00000E+01 


Figure 8 


31 









SESSION 53 - WORK-IN-PROGRESS II 



N89- 24630 


Influence of Time Lag and Noncolocation on 
Integrated Structural/ Control System Designs^ 

R. A. Manning* and L. A. Schmit** 


t 


This research was supported in part by NASA Research Grant NSG 1490. 


* Member of the Technical Staff, 

TRW Space and Technology Group 

** Professor of Engineering and Applied Science, 
University of California, Los Angeles 


i 


PRECEDING PAGE BLANK NOT FILMED 


35 


Introduction 


Recent research efforts have led to the development of simultaneous structural/control system 
design procedures [1 — 3]. Absent in any of this work is the time delay present in the control system 
sensors and actuators and the computational time delay for synthesizing actuator commands from 
sensor measurements. Madden [4] has shown that the time delay present in the control system can 
have profound effects on the resulting system performance and stability regardless of its source. In 
addition, many of the simultaneous structural/control system design procedures have used colocated 
sensors and actuators for implementation of the control system. In actual practice, colocation is not 
always possible (e.g., actuator output forces are based on optical quality measurements such as line- 
of-sight). Spector and Flashner [5] and Bong Wei [6] have raised the issue of stability degradation 
when using noncolocated sensor and actuators. 

This work extends the integrated structural/control system design procedure reported in Refer- 
ence [3] to include the effects of time lag and noncolocation of sensors and actuators on the resulting 
optimum designs. 

Optimum Design Problem Statement 


In Reference [3] the integrated controls/structure optimum design problem was posed as either 
the mass minimization problem given in equations (1) through (4) or the control effort minimization 
problem given in equations (5) through (8). The set of behavior constraints, g m (d,t), consists of 
time parametric upper bounds on the peak transient dynamic displacements and accelerations 
at selected degrees-of-freedom as well as upper or lower bounds on selected natural frequencies. 
The vector of design variables, d , consists of finite element box beam cross-sectional dimensions, 
spherical nonstructural mass element radii, and nonlinear on/off control system velocity thresholds 
and actuator output force magnitudes. 

The solution to either the mass minimization problem or the control effort minimization problem 
is found by solving a sequence of explicit approximate problems. Each approximate problem is 
constructed using first order hybrid approximations for all of the critical (or near critical) behavior 
constraints as well as for the control effort objective function in (5) or upper bound constraint in (3). 
Time parametric peak transient dynamic displacement and/or acceleration sensitivities for use in 
the hybrid approximations are calculated in an efficient manner using the Wilkie-Perkins essential 
parameter sensitivity method [7, 8] which significantly reduces the amount of time stepping needed 
to obtain these sensitivities. 


min W(d) 

(1) 

subject to g m (d, t) > 0 

(2) 

E{d,t f ) < E u 

(3) 

d/ < dj < dj U 

(4) 

min E(d, tf) 

(5) 

subject to g m (d,t ) > 0 

(6) 

W(d) < w u 

(7) 

4j <dj< i,‘ 

(8) 
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Numerical Example 


The 21 degree-of- freedom aluminum grillage structure shown is used to examine the effects of 
control system delay and noncolocation of sensors and actuators on the resulting optimum designs. 
Nine analysis box beam finite elements are linked to yield five design elements. The external load 
shown was applied at node 8 of the structure so as to excite both cantilever bending modes and 
torsional modes. A colocated sensor/actuator pair is located at node 6 of the structure to try 
and reduce dynamic response. Upper bounds of 9.0 X 10~ 4 m were placed on the peak dynamic 
displacement response at nodes 5, 6, and 7. All dynamic response calculations were carried out for 
1 second using 10 retained modes (frequency content up to 100 Hz), 2% modal damping, and a time 
step of 0.005 seconds. Both the minimum weight problem (with an upper bound placed on control 
effort) and the minimum control effort problem (with a weight cap) were used to demonstrate the 
effects of control system time lag on the resulting optimum designs. Results for this problem were 
reported in Reference [3] using colocated sensors and actuators and no time lag. 


O node 



P 8 (t) (N) 
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Numerical Results - With Time Lag 


The actuator output forces were represented as exponential growth functions with the value of 
the time lag defined as the time it takes to generate 90% of the maximum actuator output force 
from the time the actuator is commanded to generate force. For each value of time lag, a complete 
optimization was done from the same starting point. The results of the optimizations for each time 
lag are shown in the accompanying tables and figures. 

The figures show that the value of the objective function increases with increasing time lag when 
displacement constraints drive the optimum design. Also plotted in the figures are the corresponding 
values of the velocity threshold design variable, e;. One can see that as the time lag increases, the 
value of the velocity threshold at the optimum design decreases. This indicates that the optimizer 
is compensating for the time lag by commanding the actuators to generate force sooner. At time lag 
values greater than 15 milliseconds the minimum weight increases nearly linearly since the velocity 
threshold design variable is against its lower bound side constraint. 

Time histories for the critical constraints and corresponding actuator forces are shown on the 
following pages. 


Table 1: Effect o: 

■ Time Lag on the Optimum Designs 

Case 

Objective 

Function 

Time Lag 
(sec) 

Objective 

Function 

Value 

Critical 

Behavior 

Constraints 

1A 

mass 

0.0001 

346 


IB 

mass 

0.0010 

357 

Xy7 

1C 

mass 

0.0020 

381 

Xy7 

ID 

mass 

0.0030 

404 

Xy 7 

2A 

effort 

0.0001 

12.75 

X y6 ,X y7 ,W 

2B 

effort 

0.0010 

14.65 

XyS , Xy7 , W 

2C 

effort 

0.0020 

15.85 

Xye,X y7 ,W 

2D 

effort 

0.0030 

16.15 

Xy 6 ,X y7 ,W 
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Numerical Results - With Time Lag 
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CASE 1 A 


CASE 1C 









Numerical Results - Noncolocated Sensors and Actuators 


The minimum weight design problem was run again using noncolocated sensors and actuators 
without time lag. In this case, simple estimators (based on kinematics) had to be constructed to 
generate velocities at actuators locations from the measured velocities at sensor locations. Case 3A 
in the table below gives the result obtained in Reference [3] for colocated sensors and actuators. 
The two configurations of noncolocated sensor/ actuator locations, their respective estimators, and 
the resulting optimum designs are given in the table as cases 3B and 3C. 

Case 3B and 3C time histories for the critical constraint at the optimum designs, X y5 , along with 
the actuator force outputs are shown in the accompanying figures. In both instances, peak dynamic 
displacements are reduced in the time period considered when compared with the uncontrolled 
case. However, towards the end of the time history, both cases show high frequency oscillations 
superimposed on the response. In fact, both responses are unstable if carried out for a longer period 
of time. This instability is caused entirely by the fifth mode (17.8 Hz in case 3B, 18.2 Hz in case 3C) 
where the sign of the translational component of the eigenvector at the sensor location is opposite 
to the sign at the actuator locations. As pointed out by Bong Wie [6], systems with noncolocated 
sensors and actuators tend to be closer to the stability bounds for just that reason. The actuator 
force time histories show that the control system is pumping energy into the system at exactly the 
frequency of the fifth mode, thus driving it unstable. 


Table 2: Effect of Noncolocation on the Optimum Designs 


Case 

Nodal 

Sensor 

Locations 

Nodal 

Actuator 

Locations 

Estimator 

Equations 

Objective 

Function 

Values 

3A 

5,7 

5,7 


259 

3B 

6 

5,7 

Xy$ ^ Xy£ + 5 X$ x $ 

Xyj X y 6 — 5Xe x e 

406 

3C 

3 

5,7 

XyS = Xy3 + 5X^3 + 5X$ z3 
X y 7 = X y3 — 5 Xo x3 + 5X<9 z3 

423 
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Numerical Results - Noncolocated Sensors and Actuators 


CASE 3B 


CASE 3C 











Concluding Remarks 


The effects of including time lag and allowing noncolocated sensors and actuators on optimum 
designs has been explored. Results of this study show that neglect of time lag in the control system 
can lead to unconservative designs (i.e., lower objective function designs than can be physically 
realized). However, the time lag results indicate that it is feasible to incorporate this refinement 
within a design optimization procedure which includes direct constraints on peak transient dynamic 
displacements at specified degrees of freedom. Furthermore, it is found that the optimization 
procedure compensates for the presence of time lag in the system by lowering the velocity thresholds, 
thus turning on the control system sooner. On the other hand, the use of noncolocated sensors and 
actuators can result in convergence to a dynamically unstable system when trying to control modes 
where the sensors and actuators are out of phase with respect to one another. The results for 
the noncolocated sensors and actuators indicate that it will be necessary to add constraints on 
appropriate dynamic stability measures in order to prevent unstable behavior. 
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Summary 


A structural optimization algorithm was researched including 
global displacements as decision variables. The algorithm was applied 
to planar reinforced concrete frames with nonlinear material behavior 
submitted to static loading. The flexural performance of the elements 
was evaluated as a function of the actual stress-strain diagrams of 
the materials. Formation of rotational hinges with strain hardening 
were allowed and the equilibrium constraints were updated accordingly. 
The adequacy of the frames was guaranteed by imposing as constraints 
required reliability indices for the members, maximum global 
displacements for the structure and a maximum system probability of 
failure. 


Previous Research 


Structural frame optimization problems have been usually 
formulated based on the cycling between two distinct phases: analysis 
and optimal design. The option described in this work combines both 
phases by the addition of the global displacements to the set of 
design variables, option researched by several authors (ref. 1 and 2). 
The main purpose of this strategy was to determine the benefits of 
extending the linear static formulation to nonlinear static structural 
problems using the secant stiffness method. The reason behind this 
research was the fact that the global stiffness changes created by the 
nonlinear behavior would be considered simultaneously with the changes 
of the element sizes, thus improving convergence. 

The first step of the research was to optimize elastic plane 
frames with elements with rectangular sections submitted to static 
loading. The objective function was the volume of the structure and 
constraints of the optimization problem were equalities representing 
global equilibrium of the structure and inequalities for the limits on 
global displacements and the maximum flexural stresses. The strategy 
adopted consisted of transforming the constrained problem in an 
unconstrained one using the method of the Augmented Lagrangian 
Multipliers (ref. 3). The unconstrained minimization was solved using 
the Hooke and Jeeves method. 

The results with this formulation were encouraging and the 
optimal solutions were found. The convergence rate was dependent on 
the initial design, scaling, penalty parameters and lagrangian 
multipliers values. The computational effort was considerable when 
compared with other explored techniques based on optimality criteria 
and mathematical programming methods. To improve the efficiency of 
the algorithm a gradient technique was implemented to solve the 
unconstrained minimization problem. This improvement was unsuccessful 
since the Augmented Lagrangian function was very steep, with large 
sensitivity to any small variation of the displacement variables. 
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Nonlinear Formulation 


The following logical objective was to extend this strategy to 
nonlinear reinforced concrete frames. The typical frame element has 
rectangular cross section and is doubly reinforced with equal amount 
of flexural steel on both sides. The model adopted for the inelastic 
reinforced concrete element was the one component model , where 
rotational springs are added to the ends of the elastic element to 
simulate the formation of plastic binges at the extremities of the 
element (ref. 4). The stiffnesses of the linear elastic element 
stiffness and the springs was condensed using the flexibility 
formulation. 

The determination of the characteristics of each reinforced 
concrete section was based on the stress strain diagrams for the 
concrete and the reinforcing steel . The yielding and ultimate moments 
for each cross section were used to determine the characteristics of 
the springs for each element. The spring stiffness was considered 
infinite whenever the element moment was below the yielding moment. 
When the moment was above the yielding value the spring stiffness was 
updated and, since there was no incremental loading or unbalanced 
iteration of the structure, the secant stiffness was adopted for the 
spring stiffness (ref. 5). 


ONE COMPONENT MODEL AND MATERIAL CHARACTERISTICS 


Linear Elastic Eleaent 



One Component Reinforced Concrete Element 
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Procedure Implementation 


The objective function was the cost of the materials: concrete 
and steel. The constraints included the equilibrium equalities, 
maximum global displacements and maximum probability of failure for 
each element. The equilibrium constraints were evaluated every time a 
design variable changed with the corresponding updating of the spring 
stiffnesses. The values of the maximum displacements were dictated by 
serviceability constraints like maximum joint rotations or story 
drifts. The maximum element probabilities of failure were chosen with 
current practices of structural design codes (ref. 6). 

Flexural element actions are the most important in small and 
medium sized frames for the definition of section sizes and 
longitudinal steel. The evaluation of element reliability was based 
on the corresponding flexural failure function (ref. 7). The same 
approach was used to evaluate the system probability of failure. The 
basic variables considered were the compressive strength of concrete 
and the external loads. 

The system probability of failure was evaluated at the 
mechanism level at the end of each optimization cycle. If the value 
of the system probability of failure was not satisfactory the 
optimization was restarted using a different limit of element 
probability of failure for the elements involved in the failure 
mechanism. The system probability of failure was obtained using the 
beta-unzipping method (ref. 8). In summary, the elementary mechanisms 
of failure were determined using Watwood's method (ref. 9) and the 
correspondent failure functions formed. These mechanisms were then 
combined linearly and the related probabilities of failure calculated, 
while rejecting those combinations with values outside given 
intervals . 


SPRING SECANT STIFFNESS 


^ Honent 

Spring Monent-Rotation Oiagran 
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Hy - Yielding aoient ry - Yielding rotation 
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Optimization Results 


The first approach to solve the optimization problem used the 
Augmented Lagrangian function and the Hooke and Jeeves method as the 
unconstrained minimization technique. Considerable effort was put 
into this formulation with several options for the starting points, 
combinations of penalty parameters, scaling techniques and number of 
cycles. Three structures were tested with different levels of 
complexity. 

In some cases the values obtained were close with those 
corresponding to the expected optimal values. Reliability constraints 
were satisfied, displacements were within the limits and equality 
constraints were satisfied. However, convergence was difficult to 
obtain and largely dependent on several different choices made at the 
start of the optimization cycle. At the same time the element forces 
were not in accordance with the assumed secant spring stiffnesses 
showing lack of convergence of the nonlinear iteration process. 

To improve convergence of the nonlinear equilibrium of the 
structure an intermediate phase was created in the optimization 
process. The displacements were removed from the optimization cycle. 
This phase corresponded to the solution of the equilibrium equations 
every time a cycle of variable optimization in the Hooke and Jeeves 
was completed. The displacements were obtained from the equilibrium 
equations assuming the values of the secant spring stiffnesses as 
those at the end of the cycling optimization. The results were 
nevertheless the same as before and for that reason another technique 
was implemented. 


EXAMPLES TESTED 
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The Generalized Reduced Gradient method was chosen because of 
it's characteristic of solving iteratively a set of nonlinear 
equations. The algorithm used (ref. 10) performed very well for the 
elastic case with convergence in most of the cases. It proved to be 
almost insensitive to the initial design points. The extension to the 
nonlinear material behavior is however in process. A first phase of 
this extension consisted of assuming for the secant spring stiffness a 
yielding value, i.e., corresponding to the ratio of the yielding 
moment and the yielding rotation whenever the element moment was 
greater than the yielding moment. Equilibrium constraints were 
satisfied, element moments were in accordance with the assumed spring 
stiffnesses and the global displacements were in accordance with 
assumed spring stiffness values. 

The second phase of transforming the spring stiffness from the 
yielding stiffness to the secant stiffness is presently being 
researched. The initial design for the values presented using the 
secant spring stiffness formulation was obtained with the elastic 
stiffness version having as ultimate moment the yielding moment. The 
secant stiffness values were limited to a minimum value corresponding 
to the ratio of the ultimate moment and the ultimate rotation to 
prevent severe oscillations of these values. The results verified the 
equilibrium constraints within certain tolerance, the global 
displacements were those corresponding to the element spring 
stiffnesses, reliability constraints were satisfied and there was an 
improvement of the cost of the structure. All elements have yielded 
corresponding to the expected results from a optimal configuration. 
However, there are discrepancies in the moments at the joints which 
shows a insufficient convergence of the equilibrium constraints. 


RESULTS WITH YIELDING STIFFNESS AND SECANT STIFFNESS 


100 ’ 



O - Spring 



ELE 

MENT 

INITIAL 

ELASTIC 

YIELDING 

SECANT 

b 

h 

As 

b 

h 

As 

b 

h 

As 

b 

h 

As 

1 

5.0 

10.0 

1.0 

2.0* 

6.0* 

.20 

2.0* 

6.0* 

.21 

2.0* 

6.0* 

.38 

2 

5.0 

10.0 

1.0 

2.0* 

6.0* 

.23 

2.0* 

6.0* 

.21 

2.0* 

6.0* 

.25 

3 

5.0 

10.0 

1.0 

2.0* 

11.1 

.73 

2.0* 

9.75 

.64 

2.0* 

9.81 

.64 

4 

5.0 

10.0 

1.0 

2.0* 

12.4 

.B1 

2.0* 

10.9 

.71 

2.0* 

10.9 

.72 

Cost 

18,000 

6,890 

6,296 

6.505 


b - base (in); h - height (in); As - steel area (in2); * - lower bounds. 
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Directions Headed 


The Augmented Lagrangian function together with the Hooke and 
Jeeves for the unconstrained minimization although insensitive to 
discontinuities of the constraints is not the best choice when 
compared with the Generalized Reduced Gradient method. More research 
has to be done on the robustness of the convergence of the nonlinear 
iteration process (ref. 11). 

A possible improvement would consist of a similar enhancement 
to the one applied in the Augmented Lagrangian formulation with an 
intermediate solution of the displacements during the optimization 
cycle. Another possible improvement would be the use of a cycling 
procedure where the spring stiffness values were kept constant during 
the optimization and updated at the end with consequent optimization 
cycle until there was stabilization of the spring stiffness values. 

The integrated approach proved itself adequate for the elastic 
stiffness if an adequate mathematical programming technique is chosen. 
It is logical to expect that it will probably perform well in the case 
of inelastic stiffness if adequately integrated with some kind of 
nonlinear structural analysis technique. 
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INTRODUCTION 


Sensitivity analysis is fundamental to the solution of structural optimization problems. 
Consequently, much research has focused on the efficient computation of static displacement 
derivatives (Ref. 1). As originally developed, these methods relied on analytical 
representations for the derivatives of the structural stiffness matrix (K) with respect to the 
design variables ( bj ). To extend these methods for use with complex finite element 
formulations and facilitate their implementation into structural optimization programs (eg. Ref. 
2) using general finite element analysis codes (Refs. 3-4), the semi-analytic method (Refs. 5-6) 
was developed. In this method we approximate the matrix <9K/ dbj by finite difference. 

Although it is well known that the accuracy of the semi-analytic method is dependent on 
the finite difference parameter, recent work (Ref. 7) has suggested that more fundamental 
inaccuracies exist in the method when used for shape optimization. Another study (Ref. 8) has 
argued qualitatively (for the case of a cantilevered beam) that these errors are related to non- 
uniform errors in the stiffness matrix derivatives. 

In the following we will investigate the accuracy of the semi-analytic method. We first 
develop a general framework for the error analysis and then show analytically that the errors in 
the method are entirely accounted for by errors in AK/ A b { . Furthermore, we demonstrate that 
acceptable accuracy in the derivatives can be obtained through careful selection of the finite 
difference parameter. 


Static displacement derivatives: 
a57~- k i^ 7 u — k p ’ 

AK 

In the semi-analytic method we approximate by finite differences 

to compute the approximation = -4^ efficiently 


This gives the semi-analytic formula 


Au 
A b: 


- -K 


-i AK 


u = -K -1 p . 


Ab ■ 


which has 


been used successfully for sizing optimization. 


Barthelemy & Haftka - demonstrated large errors for shape optimization 
Pedersen, Cheng, & Rasmussen - some analysis of these errors 


Figure 1 
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ACCURACY OF THE SEMI- ANALYTIC METHOD 


To characterize the errors associated with the semi-analytic method let us examine the 
expression for Au/ A bj in terms of the approximate pseudo-load vector pS (see figure 2) and 

consider two cases. In both cases K is separated into two parts (K and K(bj )) which are 

independent of and dependent on the design variable bj , respectively. In the first case we can 
factor K(bj ) into a constant matrix K bj and a scalar function f(bj ). As a result, the approximate 

pseudo-load vector is a simple scaling of the true pseudo-load vector pS and the semi- 

analytic method yields displacement derivatives which are scaled with respect to the analytic 
derivatives. In this case the accuracy of the derivative is only dependent on the accuracy of 
Af/A bj . 

If, as is often the case, K(bj ) can not be factored as described above then the errors in the 
displacement derivatives may have a significantly different form. In this case the approximate 
pseudo-load vector is not a simple scaling of the true pseudo-load vector. Geometrically, this 
means that both the shape and length of the approximate pseudo-load vector may be incorrect. 
Also, since pS is a function of both the error matrix Ej and the displacement field, u, the 

accuracy of the derivatives may depend on the number of elements in the structural model and 
the location, within the model, of the element(s) dependent on bj . 


Derivative given by = 


• K factorable in b { : K = K + fib ,.)K 


b. 


* A / (b ) 

Pi = ~Kb7~ K b. u = 0 + £,• s >~db~^b . u = ( 1+ e X- 


Thus the derivatives scale: 


Au 

A b ; 


= ( 1 + e. 


3D 


• K noj factorable in b ( : K= K + K(b t ) 


s AK (b .) 

P,= 


A b; 


u = 


dK(b.) 


db ; 


— + E . 


u = p* . + E .u 

r i i 


Derivatives do not scale: ^-= -K ^E.u) | 


Figure 2 
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EXAMPLE - CANTILEVERED BEAM 


To illustrate these ideas, consider a cantilevered beam modeled as an assemblage of 
beam type finite elements. The element level stiffness matrix for the n-th element is shown in 
figure 3. Clearly, K e (and thus K) is factorable for the element height Qi) and width ( zv ) 
variables but is not factorable for the length variable (/). Quantitatively then (based on our 
previous arguments) we expect the following when using the semi-analytic method: since w 
appears linearly, the derivatives of the displacements with respect to w will be exact and the 
relative error in the derivatives of the displacements with respect to h will be uniform and 

depend only on the accuracy of A/j 3 /A/i. However, the relative error in the displacement 
derivatives with respect to / may be non-uniform and may depend on the number of elements 

used to model the beam as well as on the accuracy of AK/ A/. To confirm this, we will now 
derive analytical expressions for the relative error in these displacement derivatives. 


Consider the finite element formulation for a beam element: 


K 


e _ 
n — 



' 6 

~lln 

-6 

-3 l n 

2 EI n 

: 

it 

3 /„ 

l 2 

1 n 

ll 

sym 

: 

6 

3 In 

2 1 2 
^ l n 


where: 




12 


Clearly, K c rt may be factored for the sizing variables h and w but not 
for the shape variable / 


Figure 3 
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ERROR ANALYSIS FOR THE CANTILEVERED BEAM 


Consider the initially uniform cantilevered beam of length L shown in figure 4. The 
beam has a rectangular cross section and is subject to concentrated force (F) and moment (M) 
loadings at the tip. In order to investigate the accuracy of the tip displacement derivatives we 
can derive analytical expressions for the semi-analytic derivatives and compare them to the 
known true derivatives (Ref. 9). We begin with the expressions for the displacements and 
rotations along the length of the beam as shown below. Now let the beam be composed of N 
elements of length l=L/N numbered from 1 to N, starting at the root. If the nodes are numbered 
from 0 to N, starting at the root, then the n-th node is located an x=nl. Substituting for x in the 

equations for u and 6 yields a set of discretized equations for u and 6. To complete the 
derivation we need expressions for the entries of associated with the tip displacement d.o.f. 
and for the derivatives of the stiffness matrix with respect to the design variables b n . The 
stiffness matrix derivatives are easily derived from the expression for the element level stiffness 
matrix shown previously. The necessary entries of K' 1 can be obtained by differentiating the 
displacement vector with respect to the applied force (F). Substitution of these expressions into 
the equation for the tip displacement derivative will yield the desired analytical expression for 
the semi-analytic derivative. 


We can determine analytical expressions for the errors introduced by finite 
differencing in the S-A method by using the exact beam element formulation. 



Discretizing for a beam of N elements we can compute: 


AtF 

Ab n 


= -{k- 1 } AK 


A b„ 


u 


where {k -1 } 


t 


du 

dF 


Figure 4 
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ERROR ANALYSIS FOR THE CANTILEVERED BEAM 


In figure 5 the expressions for the semi-analytic tip displacement derivatives with 
respect to the element heights and lengths (in terms of the exact derivatives) are shown. The 
expressions have been simplified to the case where F= 0. For h, the stiffness matrix is factorable 

and, as was predicted, the relative error (e) depends only on the finite difference parameter (c) 
and is the same for all elements making up the beam. For the element lengths, K is not 
factorable and the relative error is non-uniform. In this case the relative error depends not only 
on c, but also on the element number (n) and the number of elements (N) used to model the 
beam. Increasing either n or N will cause the relative error to become larger, while decreasing 
the value of c will give better accuracy. In Refs. 7 and 8 the tip displacement derivatives with 
respect to L (Aid/AL) are investigated. This quantity is based on perturbations of all elements 
in the beam such that the quantity AL is distributed evenly among all elements. In this case the 
relative error in A u 1 / AL is equivalent to the average error (e aV g) in Aid/A/ n . For small values of 

c, E aV g is approximately proportional to cN 2 . Note that, as would be expected, in all cases the 
relative error approaches zero as c approaches zero. 


• Factorable sizing variable h : k= K+/(/i 3 )K a 



• Non-factorable shape variable l : K = K + K(/„) 



Figure 5 
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COMPUTATIONAL RESULTS: SHAPE VARIABLES 


To demonstrate the analytical results presented previously the cantilevered beam was 
modeled with beam elements and the tip displacement derivatives with respect to the element 
lengths were calculated numerically via the semi-analytic method. In the figure below these 
derivatives (for c=.01) are plotted (normalized by the true derivative) versus the element 
number for four different beam discretizations (5,10,15 and 20 elements). The numerically 
generated data points are represented by the symbols shown on the plot. The analytic results 
appear as the underlying curves. Note that the computed values are in complete agreement 
with the analytic values. As predicted, the relative error depends on both n and N. For this 
problem, the error increases with N and increases as we move along the beam from the root to 
the tip. 

In figure 6 we also show the average error of the tip displacement derivatives, with 
respect to the element length, as a function of the number of elements in the beam. In this case 
the computed and predicted values are plotted for various values of c. In addition, the 
equivalent numerical data from Ref. 7 (represented by the square symbols) is also shown. 
Note, again, the excellent agreement between the computed and predicted errors. Clearly, the 
error decreases rapidly as c is decreased. For c=.00001 the average error is less than 1.0%. for 
N= 20. 


Average Error of Tip 

Tip Displacement Derivatives With Displacement Derivatives 

Respect To Dement Length (c=.001) With Respect to Element Length 




Figure 6 
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ERROR ANALYSIS FOR NON-PLANAR BEAM 


We have seen that for the cantilevered beam problem uniform errors occur for the tip 
displacement derivatives with respect to sizing design variables and nonuniform errors occur 
for shape variables. However, our general error analysis predicts that the nature of the error 
depends on the factorability of K and not necessarily on the classification of the design variable. 
We will now show that nonuniform derivative errors can (and usually do) occur for sizing 
design variables. 

Consider the case of a rectangular cantilevered beam where the principal axes of the 
beam elements are oriented at some angle to the global coordinate system. In this case the 
stiffness matrix is not generally factorable for w and h since both variables contribute to the 
stiffnesses in the global coordinate system. As a result nonuniform error in the global 
displacement derivatives may occur. To investigate this analytically, consider the equivalent 
system shown in figure 7 where we want to calculate the derivative of the displacement at 

some angle 0 to the beam's principal axes (Ad f /A b n ). Using our previous analysis for 
displacement derivatives in the principal coordinate system we can develop the expression 

shown below for the relative error (e M ) associated with AdVA b n (Ref.9). As shown, e n is a 
function of the relative errors associated with the principal displacement derivatives Generally, 
they combine such that e n will differ for each element. Under certain conditions, however, 
uniform errors will occur. Clearly, this will be the case when 0 is some multiple of re / 2. 
Uniform error will also occur when the loading in y and z directions are related by a scaling 
factor since the resulting displacements and displacement derivatives will also be simply 
scaled. 



The relative error of the displacement derivative is given by : 


£» = 


Ad* 
A b n 


dd 

db„ 


e“-^-cos0 + d'-lr-sin 9 


db. 


db n 


dd 

db n 


cos 9 + sin 9 


• Uniform error when 9 = in/2 or P z = yP^, 

• Non-uniform error otherwise 


Figure 7 
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NONUNIFORM ERRORS - SIZING VARIABLES 


To numerically illustrate that nonuniform errors can occur for displacement derivatives 
with respect to sizing design variables we calculated the derivative of the tip displacement for a 
cantilevered beam with respect to the element heights using the semi-analytic method. The 
beam was rotated so that its principal axes were oriented at 45 degrees relative to the global 
coordinate system. A tip force and moment were applied parallel to the global axes. In figure 
8 the derivatives are plotted (normalized by the analytic derivatives) as a function of the 
element number for c=0.01 and various beam discretizations (5,10,15 and 20 elements). As 
expected, the derivatives depend on the element number and the number of element used to 
model the beam. In this case the relative error decreases with increasing N and decreases as we 
move along the beam from the root to the tip. 

In figure 8 we also show the normalized tip displacement derivatives plotted versus the 
element number for N= 20 and various values of c. Note that the relative error decreases 
rapidly as c is decreased. Also, the magnitudes of the errors are significantly less, for a given 
value of c, than we found for the derivatives with respect to the element lengths. For N= 20, 
acceptable accuracy is obtained for values of c as large as 0.01. 


Tip Displacement Derivatives With 
Respect To Element Height (c=.0l) 


Tip Displacement Derivatives With 
Respect To Element Height (N^20) 




Figure 8 
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ERROR ANALYSIS - AUTOMOTIVE FRAME 


To study the accuracy of the semi-analytic method for more practical problems, consider 
the half model of an idealized automobile frame structure shown in figure 9. The structural 
model consists of 33 three-dimensional beam-type finite elements each having a rectangular 
cross section. The structure is simply supported at the front suspension attachment points (A) 
and loaded in the vertical direction at the rear suspension attachment points (B). Boundary 
conditions are applied to the center line grid points to enforce an anti-symmetric structural 
response. The net effect of the loading and boundary conditions is to cause torsion of the 
structure about the centerline. In this case we calculated the semi-analytic derivatives of the 
vertical displacement at point C with respect to the thickness, width, height and length of each 
element in the structure and compared them against the analytic derivatives. In all cases the 
accuracy of the derivatives varies from element to element. The results of the comparison are 
summarized in the plot below. For each type of design variable (length, height, width and 
thickness) the minimum, maximum and average errors in semi-analytic derivatives are plotted 
as a function of the finite difference parameter (c). Note that each type of design variable 
exhibits a different level of accuracy, for a given value of c, with the length variable being the 
worst and thickness being the best. This can be attributed to the varying degrees of 
nonlinearity of the stiffness matrix with respect to these variables. For a thin walled box beam 
the section properties are nearly linear functions of t and therefore the accuracy of the 
displacement derivatives with respect to t is much better than that for b, h and /. In general, the 
careful selection of design variables or other intermediate variables (e.g., beam section 
properties) for the derivative calculations will yield more accurate derivatives for any given 
value of c. 



RELATIVE ERRORS AS A FUNCTION OF FINITE DIFFERENCE 
PARAMETER AND DESIGN VARIABLE TYPE 



Figure 9 
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SUMMARY 


The inaccuracy of the semi-analytic method for computing static displacement 
derivatives for both shape and sizing design variables has been shown to be the result of errors 
in the pseudo load vectors. Two types of errors were identified. In the first case the errors in 
the finite difference approximation to the stiffness matrix derivatives resulted in a scaling of 
the pseudo load vector which, in turn, causes the derivatives to be uniformly scaled relative to 
their true values. In this case the magnitude of the error depends only on the finite difference 
parameter, c. In the second case, errors in the finite difference operation lead to a distortion of 
the pseudo load vectors and nonuniform errors in the displacement derivatives. These errors 
may be dependent on the location (within the structure) of the element(s) associated with the 
design variable and the discretization of the structure, as well as c. 

The results of the error analysis were demonstrated numerically for a cantilevered beam 
and an idealized automobile frame structure. It was observed that for a given value of c that 
the errors in the derivatives for shape design variables were significantly larger than those for 
sizing variables. However, in both cases the relative errors could be adequately controlled 
through the proper choice of the finite difference parameter. It should be noted that relatively 
small values for c may be required to compute sufficiently accurate derivatives. This suggests 
that it may be necessary to compute the finite difference approximations to the stiffness matrix 
derivatives in double precision to avoid roundoff errors. Also, by carefully choosing 
intermediate variables which appear linearly (or nearly so) in the stiffness matrix, greater 
accuracy in the the finite difference approximation can be obtained. 


• Errors reported for the semi-analytic method have been shown to 
be due to errors in the finite difference approximation of the 
stiffness matrix derivatives 

• We can adequately control errors by careful choice of the 
finite difference parameter 

• Errors may occur for both shape and sizing variables 

• For a given value of the finite difference parameter, errors in 
the derivatives of the shape variables were larger than the 
sizing variables 


Figure 10 
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ABSTRACT 


An analytical method is being developed to determine the signature of an 
acoustic emission waveform from a growing crack and the results of this analysis 
are compared to experimentally obtained values. Within the assumptions of linear 
elastic fracture mechanics, a two dimensional model is developed to examine a 
semi- inf inite crack that, after propagating with a constant velocity, suddenly 
stops. The analytical model employs an integral equation method for the analysis 
of problems of dynamic fracture mechanics. The experimental procedure uses an 
interferometric apparatus that makes very localized absolute measurements with 
very high fidelity and without acoustically loading the specimen. 

INTRODUCTION 


Acoustic emission testing is a method of nondestructive evaluation that 
detects stress wave emissions from fracture and deformation processes within a 
loaded body. This testing differs from other methods of nondestructive 
evaluation in that the signal being detected is released from within the specimen 
rather than being created by the nondestructive testing method. The technique 
offers a distinct advantage over more conventional nondestructive testing 
techniques because it allows for the real time monitoring of in-service 
structures. Some of the potential source mechanisms of acoustic emission 
include: crack propagation and arrest, fretting among fracture surfaces, 
dislocation movement, microcracking, twinning and phase transformations. In 
addition to these failure related mechanisms, other phenomena such as fastener 
fretting, structural vibration and electromagnetic noise can create spurious 
signals which are detected by the acoustic emission instrumentation. Of 
fundamental importance for the advancement of the current state of acoustic 
emission technology is the isolation and identification of the signal from a 
growing crack. The technology for detecting and locating internal sources of 
acoustic emission is well established. However, acoustic emission signals 
contain a vast amount of additional information about the source of the emission 
and the condition of the material being examined. The signal is not only 
influenced by its source but also by the specimen geometry (which effects the 
stress wave propagation from the source to the sensor) and the charac teristics 
of the sensor. A thorough understanding of each of these factors is necessary 
in order to accurately interpret the acoustic emission signature. The proposed 
solution procedure will attempt to apply methods from dynamic fracture mechanics 
and wave propagation to the quantitative characterization of acoustic emission 
signals . 

This work complements previous studies by providing a development of the 
analytical form of an acoustic emission waveform caused by a crack growth event. 
An advantage of the proposed analysis is that the source for the acoustic 
emission signature is an actual crack propagation event and not a simple point 
source model. The propagation of the crack greatly influences the stress field 
in the vicinity of the crack tip, causing stress wave fronts to radiate into the 
body and on the crack surface. Acoustic emission testing detects these stress 
waves at the body's surface and relates the signal back to the corresponding 
crack propagation event. The proposed method uses an integral equation technique 
for the analysis of problems of dynamic fracture mechanics developed by Jacobs 
and Bieniek [1] . The problems of dynamic crack propagation have been the subject 
of numerous investigations in the past several years. A majority of the work 


68 



has been summarized in review articles by Achenbach [2, 3], Freund [4, 5], Rose 
[6] and in the book by Kanninen and Popelar [7], The preceding works have 
primarily been concerned with the determination of the dynamic stress field in 
the vicinity of the moving crack tip. They did not examine the effect of the 
propagation and arrest of the crack tips throughout the entire body. Freund [8] 
determined the pressure discontinuity radiated out from a crack tip when the 
crack, which is initially at rest, begins to grow. Rose [9] calculated explicit 
formulae for the stress discontinuities radiated by a suddenly starting two 
dimensional crack under tension for application to acoustic emission testing. 
Achenbach and Harris [10] examined the acoustic emission signals from a semi- 
infinite crack of arbitrary shape using the elastodynamic ray theory. Harris 
and Pott [11] investigated the surface motions excited by fracture processes at 
the edge of a buried crack. 

Previous investigators worked to identify the acoustic emission signal 
from a crack propagation and arrest source. Summaries of this work appears in 
Eitzen and Wadley [12], Pao [13] and in a book by the American Society for 
Nondestructive Testing [14]. Hutton, Friesel, Graham and Elsley [15] 
successfully characterized acoustic emission signals in laboratory investigations 
using statistical pattern recognition algorithms which characterize signals 
empirically on the basis of features observed in a large number of events. These 
methods are strictly empirical in nature and provide little insight into the 
fracture process. Other investigators concentrated on the geometrical effects 
of the acoustic emission signal. Pao, Gajewski and Ceranoglu [16] and Ceranoglu 
and Pao [17, 18, 19] examined the propagation of an acoustic emission signal in 
an elastic plate. These solutions are not empirical, but are based on the 
generation and propagation of elastic waves in a wave guide. The solutions, 
which examine point sources inside an infinite elastic plate, use a generalized 
ray theory and integral transform techniques. They provide numerical results 
for the surface displacements for a variety of dynamic nuclei of strains, 
including concentrated forces and couples. Individual or combinations of these 
sources are used to model the dynamic processes of material defects. Kim and 
Sachse [20, 21, 22, 23] investigated both the analytical and experimental 
signature of an acoustic emission waveform. 


ANALYTICAL METHOD 


The integral equation in the present application is in two variables, a 
spatial coordinate (x) and time (t). Within the assumptions of linear elastic 
fracture mechanics, the dynamic stresses caused by a prescribed crack growth 
event in an infinite two dimensional body are calculated. These results can be 
used to calculate displacement as a function of time at any point within the 
body. The first step of this analysis, summarized in [24, 25], uses an influence 
function to formulate an integral equation that expresses the boundary conditions 
in the plane of the crack. The steps for the calculation of the dynamic stresses 
are as follows: 

(a) Determination of the influence (or Green's) function of the problem, 
which is the dynamic displacement of an elastic half- space subjected 
to a unit concentrated impulse acting at the point of, and normal 
to, its edge. 
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(b) Formulation of the integral equation of the problem. This integral 
equation, with the influence function as the kernel and the normal 
stress in the plane of the crack as the unknown function, expresses 
the boundary conditions in the plane of the crack - a stress free 
crack surface and continuity of displacements outside of the crack. 


(c) Solution of the integral equation. 

The solution presented is for a semi- inf ini te crack that is symmetrically loaded 
(Mode I). First, solve for the influence function, U y (x-x' , t- t ' ) , in closed form 
using integral transform methods. This is accomplished by taking a one-sided 
Laplace transform in time t and a two-sided Laplace transform in x. In the 
transform space, the two uncoupled partial differential equations are replaced 
by two uncoupled ordinary differential equations. The determination of the 
inverse transformations of the required surface displacement component is 
accomplished using the Cagniard-de Hoop method. 

To formulate the integral equation, assume a crack exists at time t=0 with 
its tip located at x^a(O) and y=0 . For time t>0 , the crack moves from x=a(0) 
to x=a(t). The two relevant boundary conditions are that the newly formed crack 
faces are stress free and that the vertical displacement in front of the moving 
crack tip is zero. Both of these boundary conditions are met by: 

(a) Removing the existing known static stress, ayy^PCx) , and assuming that 
instead a new unknown time dependent stress, a yy -=F(x,t), develops. 

(b) Requiring that the new stress distribution be such that there is vertical 
displacement continuity in front of the moving crack tip. 

The continuity boundary condition can be expressed in terms of the influence 
function, U y (x-x' , t- t ' ) , as: 



P s (x')U y (T - x'J - t')dx'di' + 



*'j')U y (T-x',t-i')dx'di' = 0 


( 1 ) 


The above is a Volterra integral equation of the first kind in the variables x 
and t. To provide a simple solution of this integral equation, assume some 
spatial form of the unknown stress distribution, F(x',t'). Assume, further, 
that the spatial distribution of F(x',t') contains a square root singularity at 
its tip location a(t'), which is the same spatial form of stresses as a static 
crack with its tip located at a(t'). However, F(x',t') must contain an unknown 
time function, K(t'). Thus, the unknown stress in front of the moving crack tip 
is assumed to have the spatial form of its corresponding static crack multiplied 
by some unknown time function. It should be noted that due to the presence of 
step functions in the influence function, the infinity limits in the x' 
integration of the integral equation can be replaced by the distance that the 
fastest wave will travel in the elapsed time, t-t' . 

The direct quadrature method [26, 27] is used for the solution of the 
integral equation. There are, however, two refinements which are necessary in 
the numerical solution of the integral equation. The first refinement is a 
"subdivision" of the time step, delta t, in the evaluation of the numerical 



integrations, while the second is a mid-point product integration scheme [28] 
which is employed to handle the singularity of the kernel at t'=t. For the 
steady state case of a crack propagating with a constant velocity, the 
calculated value of K(t) is a constant that is only a function of the crack tip 
velocity. As the crack tip speed increases, the corresponding constant value 
of K(t) will decrease. The results of the case for a crack that suddenly stops 
after propagating is that the calculated value of K(t) discontinuous ly jumps to 
the value of the corresponding static stress; there is no transition zone and 
the stress never increases above the value for an equivalent static crack. 


The displacement at any point within the infinite body is determined using 
the previously calculated dynamic stress, F(x,t) and two new influence functions 
UyyCx-x' ,y , t- t ' ) and Uyy(x-x' ,y , t- t ' ) . These influence functions represent the 
horizontal displacement u* and the vertical displacement Uy, respectively, at 
point (x,y) within an elastic half-space that is subjected to a vertical unit 
impulse surface loading at x'. The solution for these new influence functions 
is accomplished using integral transform techniques and the inversion is again 
performed using the Cagniard-de Hoop method. Convolution integrals are 
developed for the vertical and horizontal displacements for any point, (x,y) , 
within the body by determining the displacements due to the application of the 
previously calculated dynamic stress distribution and the removal of the initial 
static stress distribution. The displacements are given by: 
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These integrals can be evaluated numerically. Difficulties arise in the x' 
integration due to both the integrands' complexity and the presence of 
singularities . There are singularities of different strengths and 
discontinuities associated with the moving crack tip, the original crack tip 
and the various wave fronts. Each of these singularities must be investigated 
separately to properly evaluate the integrals. To avoid numerical problems, 
the final integration technique will involve separating the x' integration 
interval into singular and non-singular regions. To evaluate the integral in 
the area of the singularity, a hybrid method is proposed. Following Davis and 
Rabinowitz in [71], the singularity is dealt with by breaking up the original 
integral in the singular region into two new integrals. One of the new 
integrals contains the singularity, but it can be evaluated analytically; the 
second integral, which is evaluated numerically, is non-singular since its 
integrand will approach zero as the potential singular point is approached. 


EXPERIMENTAL RESULTS 


The second step is to compare the results of the analytical model with 
experimentally obtained waveforms. It should be noted that the analytical 
procedure being developed is for the time prior to the arrival of stress waves 


71 


reflected from the test specimen's boundary so it is invalid for the time period 
after the fastest reflected waves interfere with the unadulterated signal from 
the crack propagation event. The experimental test procedure examines an 
existing crack that is forced to propagate a short distance and be arrested. 
The specimen used is a screw loaded, wedge opening load sample, where the crack 
unloads as it extends and the propagation is arrested before complete failure 
occurs. The specimens are made of a brittle material with a low fracture 
toughness, poly methyl methacrylate. Its optical transparency permit the size, 
geometry and location of the cracks to be readily determined. Cracks in these 
specimens are initiated by driving a sharp blade into the notch of the specimen; 
further propagation is produced by tightening the screw. 

A high sensitivity heterodyne interferometer is used to detect acoustic 
emission events. This optical device permits the high fidelity localized 
measurement of velocities from acoustic emission events arriving at various 
points on the sample surface. Since this type of measurement does not 
acoustically load the sample, the event being observed is undisturbed by the 
measurement process. The most commonly used acoustic emission sensor is the 
piezoelectric transducer. Since it must be used in direct contact with the 
specimen, the transducer will disturb the process being measured and the signal 
response will be averaged over this area of contact. Additional limitations of 
these sensors is that it is difficult to manufacture a truly broad band 
transducer, they are extremely difficult to calibrate accurately and there are 
many questions as to exactly what the transducer is measuring. 

The specimen face opposite the crack is polished and placed in the 
interferometer and becomes one mirror surface. The beam striking the face is 
approximately 1.5 mm in diameter and samples the average displacement taking 
place over this region, which is much smaller than the wavelength of the 
acoustic events being observed. The operation of the heterodyne interferometer 
is described in [30]. Briefly, single frequency laser light is split into two 
components using an acousto-optic modulator. These two components , which are 
separated in frequency by 40 MHz , are sent along two arms of an interferometer 
one of which contains the sample being monitored. The beams are recombined on 
the surface of a photodetector where they beat together at a frequency of 40 
MHz. Phase shifts in the light reflected from the sample surface result in 
equivalent phase shifts in the 40 MHz beat signal received at the photodetector. 
This carrier signal can then be demodulated to determine the time dependent 
displacement occurring at the sample surface. The detection system has a band 
width of 10 MHz which is further limited to the spectral region 0 to 2 MHz in 
order to reduce the noise in the signal. All signals are acquired on a digital 
oscilloscope and stored for later processing. The crack velocity is measured 
with conventional crack propagation gages. This also aids in determining the 
time difference between the crack growth event and the arrival of its signal at 
the measurement point. This is accomplished by pretriggering the measurement 
system on the start of the crack propagation event and not the arrival of the 
first wavefront. 


DISCUSSION 

A characteristic crack emission is shown in figure 1. Care must be taken 
to calculate the effect of wave reflections and mode conversions that occur at 
the specimen's boundaries. The experimentally obtained waveforms will be 
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interpreted using the results of the analytical model being developed in the 
first task. Since the dynamic stress calculations indicate sharp stress 
discontinuities associated with the starting and stopping phases, it is 
anticipated that there will be corresponding displacement variations that will 
become evident in the experimental modeling. Anomalies in the fracture behavior 
of the specimen included out of plane growth and some crack tunneling. The out 
of plane growth could be caused by twisting due to the bearing stress between 
the bolt and the lower crack surface. Further development of the analytical 
model is necessary before the experimentally obtained waveforms can be 
interpreted. 
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Figure 1 : Acoustic Emission Signal from a Growing Crack 
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